# cython: profile=True
cimport numpy
from numpy cimport ndarray
from numpy import zeros
from libcpp.vector cimport vector
try:
    from pulp import LpProblem, LpVariable, LpMinimize, value, LpStatus
except ImportError:
    print 'cannot import pulp'
from cython.operator cimport dereference as deref
cimport cython

cdef extern from "lpsolver.H" namespace "LPSOLVER":
    ctypedef double Real
    
    cdef cppclass Matrix:
        Matrix(size_t,size_t)
        void setval(size_t i,size_t j,Real val)
    
    cdef cppclass LpSolver:
        LpSolver(Matrix,Real *)
        void solve(Real *)


cpdef ndarray prabhu_solver_np_interface(ndarray[ndim=2,dtype=numpy.float64_t] Aeq, ndarray[ndim=1,dtype=numpy.float64_t] beq):
    cdef int nrow,ncol,i
    nrow = Aeq.shape[0]
    ncol = Aeq.shape[1]
    cdef Matrix *Aeq_c = new Matrix(nrow,ncol)
    cdef vector[double] beq_c
    beq_c.resize(nrow)
    cdef double a
    for i in range(nrow):
        for j in range(ncol):
            a = Aeq[i,j]
            Aeq_c.setval(i,j,a)
        beq_c[i] = beq[i]
    cdef LpSolver *thisptr = new LpSolver(deref(Aeq_c), &beq_c[0])
    cdef vector[double] res
    res.resize(ncol+1)
    thisptr.solve(&res[0])
    
    cdef ndarray[ndim=1,dtype=numpy.float64_t] result = zeros(res.size(),dtype='float64')
    
    for i in range(res.size()):
        result[i] = res[i]
    
    return result

cdef csolver(double** Aeq, double* beq,double* res,int nrow,int ncol):
    cdef int i
    cdef Matrix *Aeq_c = new Matrix(nrow,ncol)
    cdef double a
    for i in range(nrow):
        for j in range(ncol):
            a = Aeq[i][j]
            Aeq_c.setval(i,j,a)
    cdef LpSolver *thisptr = new LpSolver(deref(Aeq_c), &beq[0])
    thisptr.solve(&res[0])

cdef puLPsolver(double** Aeq, double* beq,double* res,int nrow,int ncol):
    prob = LpProblem('vrt',LpMinimize)
    # create variables
    w = []
    for j in range(ncol):
        newfrac = LpVariable('w'+str(j),-0.5,0.5)
        w += [newfrac]
    w += [LpVariable('wmod',0)]
    
    for i in range(nrow):
        expr1 = 0
        expr2 = 0
        for j in range(ncol):
            tmp = round(Aeq[i][j],3)*w[j]
            expr1 += tmp
            expr2 -= tmp
        prob += expr1 >= beq[i]
        prob += expr2 >= -beq[i]
    for j in range(ncol):
        prob += w[j] + w[-1] >= 0
        prob += -w[j] + w[-1] >= 0
    prob += w[-1]
    print prob
    prob.solve()
    print "Status:", LpStatus[prob.status]
    for j in range(ncol+1):
        res[j] = value(w[j])

cdef puLPsolver2(double** Aeq, double* beq,double* res,int nrow,int ncol):
    prob = LpProblem('vrt',LpMinimize)
    # create variables
    w = []
    for j in range(ncol):
        newfrac = LpVariable('w'+str(j),-0.5,0.5)
        w += [newfrac]
    
    for i in range(nrow):
        expr = 0
        for j in range(ncol):
            tmp = Aeq[i][j]*w[j]
            expr += tmp
        prob += expr == beq[i]
    print prob
    prob.solve()
    print "Status:", LpStatus[prob.status]
    for j in range(ncol+1):
        res[j] = value(w[j])

cdef puLPsolver3(double** Aeq, double* beq,double* res,int nrow,int ncol):
    prob = LpProblem('vrt',LpMinimize)
    # create variables
    f = []
    for j in range(ncol):
        newfrac = LpVariable('w'+str(j),0.,1.)
        f += [newfrac]
    
    for i in range(nrow):
        expr = 0
        for j in range(ncol):
            tmp = Aeq[i][j]*f[j]
            expr += tmp
        prob += expr == beq[i]
    print prob
    prob.solve()
    print "Status:", LpStatus[prob.status]
    for j in range(ncol+1):
        res[j] = value(f[j])
